rm(list = ls()) # Borramos variables de environment
# descomentar si es la primera vez (y requieren instalación)
# install.packages("tidyverse")
# install.packages("performance")
# install.packages("olsrr")
# install.packages("corrr")
# install.packages("corrplot")
# install.packages("skimr")
library(tidyverse)
library(rsample)
library(performance)
library(olsrr)
library(corrr)
library(corrplot)
library(skimr)Entrega I
Instrucciones (leer antes de empezar)
Modifica de la cabecera del documento
.qmdtus datos personales (nombre y DNI). IMPORTANTE: no modifiques nada más de la cabecera.Asegúrate, ANTES de seguir editando el documento, que el archivo
.qmdse renderiza correctamente y se genera el.htmlcorrespondiente en tu carpeta local de tu ordenador (pulsa el botónRendero guarda el documento conRender on Saveactivado).Los chunks (cajas de código) creados están o vacíos o incompletos, de ahí que la mayoría tengan la opción
#| eval: false. Una vez que edites lo que consideres, debes ir cambiando cada chunck a#| eval: true(o quitarlo directamente) para que se ejecuten.Recuerda que puedes ejecutar chunk a chunk con el botón play o ejecutar todos los chunk hasta uno dado (con el botón a la izquierda del anterior).
IMPORTANTE: solo se podrá subir un archivo
.htmlal campus, no se evaluarán entregas sin el.htmlcorrectamente generado. Recuerda: el código es una herramienta, no el fin en esta asignatura. Se evaluará especialmente las interpretaciones y conclusiones que detalles. Un ejercicio con el código perfecto pero sin ningún tipo de razonamiento, interpretación o conclusión no superará el 4 sobre 10 de nota.
Paquetes necesarios
Necesitaremos los siguientes paquetes (haz play en el chunk para que se carguen):
Caso práctico: análisis de datos de cáncer
El archivo de datos a usar lo cargaremos desde el csv cancer.csv.
# no cambies código
datos <- read_csv(file = "./cancer.csv")
datos# A tibble: 3,047 × 13
deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 165. 61898 260131 11.2 500. 39.3
2 161. 48127 43269 18.6 23.1 33
3 175. 49348 21026 14.6 47.6 45
4 195. 44243 75882 17.1 343. 42.8
5 144. 49955 10321 12.5 0 48.3
6 176 52313 61023 15.6 180. 45.4
7 176. 37782 41516 23.2 0 42.6
8 184. 40189 20848 17.8 0 51.7
9 190. 42579 13088 22.3 0 49.3
10 178. 60397 843954 13.1 428. 35.8
# ℹ 3,037 more rows
# ℹ 7 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>, region <chr>
Los datos representan diferentes características socioeconómicas de distintas regiones, extraídas de the American Community Survey (census.gov), clinicaltrials.gov, y cancer.gov.
¿Nuestro objetivo? Ser capaces de predecir de manera lineal y MULTIVARIANTE la variable deathRate, nuestra variable objetivo, que representa la mortalidad media de cancer, por cada 100 000 habitantes. El resto de variables son:
medianIncome: mediana de los ingresos de la región.popEst2015: población de la regiónpovertyPercent: porcentaje de población en situación de pobreza.studyPerCap: ensayos clínicos relacionados por el cáncer realizados por cada 100 000 habitantes.MedianAge: edad mediana.region: nombre de la región.AvgHouseholdSize: tamaño medio de los hogares.PercentMarried: porcentaje de habitantes casados.PctHS25_Over: porcentaje de residentes por encima de los 25 años con (máximo) título de bachillerato.PctUnemployed16_Over: porcentaje de residentes de más de 16 años en paro.PctPrivateCoverage: porcentaje de residentes con seguro de salud privado.BirthRate: tasa de natalidad.
IMPORTANTE: siempre que puedas, relaciona los valores numéricos de la salida de R con su fórmula.
Ejercicio 1 (1.25 puntos)
Chequea que no hay problemas de rango/codificación: en caso de que alguna variable tenga dichos problemas, debes sustituir los valores fueran del rango razonable por su media o mediana (sin esos datos fuera del rango), justificando por qué usas la media o por qué usas la mediana, ejecutando el código que consideres (si te atascas, sigue haciendo con el dataset original)
datos |> str()spc_tbl_ [3,047 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ deathRate : num [1:3047] 165 161 175 195 144 ...
$ medIncome : num [1:3047] 61898 48127 49348 44243 49955 ...
$ popEst2015 : num [1:3047] 260131 43269 21026 75882 10321 ...
$ povertyPercent : num [1:3047] 11.2 18.6 14.6 17.1 12.5 15.6 23.2 17.8 22.3 13.1 ...
$ studyPerCap : num [1:3047] 499.7 23.1 47.6 342.6 0 ...
$ MedianAge : num [1:3047] 39.3 33 45 42.8 48.3 45.4 42.6 51.7 49.3 35.8 ...
$ AvgHouseholdSize : num [1:3047] 2.54 2.34 2.62 2.52 2.34 2.58 2.42 2.24 2.38 2.65 ...
$ PercentMarried : num [1:3047] 52.5 44.5 54.2 52.7 57.8 50.4 54.1 52.7 55.9 50 ...
$ PctHS25_Over : num [1:3047] 23.2 26 29 31.6 33.4 30.4 29.8 31.6 32.2 28.8 ...
$ PctUnemployed16_Over: num [1:3047] 8 7.8 7 12.1 4.8 12.9 8.9 8.9 10.3 9.2 ...
$ PctPrivateCoverage : num [1:3047] 75.1 70.2 63.7 58.4 61.6 60 49.5 55.8 55.5 69.9 ...
$ BirthRate : num [1:3047] 6.12 4.33 3.73 4.6 6.8 ...
$ region : chr [1:3047] "West" "West" "West" "West" ...
- attr(*, "spec")=
.. cols(
.. deathRate = col_double(),
.. medIncome = col_double(),
.. popEst2015 = col_double(),
.. povertyPercent = col_double(),
.. studyPerCap = col_double(),
.. MedianAge = col_double(),
.. AvgHouseholdSize = col_double(),
.. PercentMarried = col_double(),
.. PctHS25_Over = col_double(),
.. PctUnemployed16_Over = col_double(),
.. PctPrivateCoverage = col_double(),
.. BirthRate = col_double(),
.. region = col_character()
.. )
- attr(*, "problems")=<externalptr>
datos |> skim()| Name | datos |
| Number of rows | 3047 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| region | 0 | 1 | 4 | 9 | 0 | 4 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| deathRate | 0 | 1 | 178.66 | 27.75 | 59.70 | 161.20 | 178.10 | 195.20 | 362.80 | ▁▇▆▁▁ |
| medIncome | 0 | 1 | 47063.28 | 12040.09 | 22640.00 | 38882.50 | 45207.00 | 52492.00 | 125635.00 | ▇▇▁▁▁ |
| popEst2015 | 0 | 1 | 102637.37 | 329059.22 | 827.00 | 11684.00 | 26643.00 | 68671.00 | 10170292.00 | ▇▁▁▁▁ |
| povertyPercent | 0 | 1 | 16.88 | 6.41 | 3.20 | 12.15 | 15.90 | 20.40 | 47.40 | ▃▇▃▁▁ |
| studyPerCap | 0 | 1 | 155.40 | 529.63 | 0.00 | 0.00 | 0.00 | 83.65 | 9762.31 | ▇▁▁▁▁ |
| MedianAge | 0 | 1 | 45.27 | 45.30 | 22.30 | 37.70 | 41.00 | 44.00 | 624.00 | ▇▁▁▁▁ |
| AvgHouseholdSize | 0 | 1 | 2.48 | 0.43 | 0.02 | 2.37 | 2.50 | 2.63 | 3.97 | ▁▁▃▇▁ |
| PercentMarried | 0 | 1 | 51.77 | 6.90 | 23.10 | 47.75 | 52.40 | 56.40 | 72.50 | ▁▂▇▇▁ |
| PctHS25_Over | 0 | 1 | 34.80 | 7.03 | 7.50 | 30.40 | 35.30 | 39.65 | 54.80 | ▁▂▇▇▁ |
| PctUnemployed16_Over | 0 | 1 | 7.85 | 3.45 | 0.40 | 5.50 | 7.60 | 9.70 | 29.40 | ▅▇▁▁▁ |
| PctPrivateCoverage | 0 | 1 | 64.35 | 10.65 | 22.30 | 57.20 | 65.10 | 72.10 | 92.30 | ▁▂▇▇▂ |
| BirthRate | 0 | 1 | 5.64 | 1.99 | 0.00 | 4.52 | 5.38 | 6.49 | 21.33 | ▂▇▁▁▁ |
#region variable cualitativa
#variable edad media se pasa
unique(datos$MedianAge) |> sort() [1] 22.3 23.2 23.3 23.5 23.9 24.2 24.4 24.6 24.7 25.0 25.2 25.9
[13] 26.1 26.2 26.3 26.4 26.5 26.6 26.7 26.8 27.0 27.2 27.4 27.5
[25] 27.6 27.8 27.9 28.0 28.1 28.3 28.4 28.5 28.6 28.7 28.8 28.9
[37] 29.1 29.2 29.3 29.4 29.5 29.7 29.8 29.9 30.0 30.1 30.2 30.3
[49] 30.4 30.5 30.6 30.7 30.8 30.9 31.0 31.1 31.2 31.3 31.4 31.5
[61] 31.6 31.7 31.8 31.9 32.0 32.1 32.2 32.3 32.4 32.5 32.6 32.7
[73] 32.8 32.9 33.0 33.1 33.2 33.3 33.4 33.5 33.6 33.7 33.8 33.9
[85] 34.0 34.1 34.2 34.3 34.4 34.5 34.6 34.7 34.8 34.9 35.0 35.1
[97] 35.2 35.3 35.4 35.5 35.6 35.7 35.8 35.9 36.0 36.1 36.2 36.3
[109] 36.4 36.5 36.6 36.7 36.8 36.9 37.0 37.1 37.2 37.3 37.4 37.5
[121] 37.6 37.7 37.8 37.9 38.0 38.1 38.2 38.3 38.4 38.5 38.6 38.7
[133] 38.8 38.9 39.0 39.1 39.2 39.3 39.4 39.5 39.6 39.7 39.8 39.9
[145] 40.0 40.1 40.2 40.3 40.4 40.5 40.6 40.7 40.8 40.9 41.0 41.1
[157] 41.2 41.3 41.4 41.5 41.6 41.7 41.8 41.9 42.0 42.1 42.2 42.3
[169] 42.4 42.5 42.6 42.7 42.8 42.9 43.0 43.1 43.2 43.3 43.4 43.5
[181] 43.6 43.7 43.8 43.9 44.0 44.1 44.2 44.3 44.4 44.5 44.6 44.7
[193] 44.8 44.9 45.0 45.1 45.2 45.3 45.4 45.5 45.6 45.7 45.8 45.9
[205] 46.0 46.1 46.2 46.3 46.4 46.5 46.6 46.7 46.8 46.9 47.0 47.1
[217] 47.2 47.3 47.4 47.5 47.6 47.7 47.8 47.9 48.0 48.1 48.2 48.3
[229] 48.4 48.5 48.6 48.7 48.8 48.9 49.0 49.1 49.2 49.3 49.4 49.5
[241] 49.6 49.7 49.8 49.9 50.0 50.1 50.2 50.3 50.4 50.5 50.6 50.7
[253] 50.8 50.9 51.0 51.1 51.2 51.3 51.4 51.5 51.7 51.8 51.9 52.0
[265] 52.1 52.2 52.3 52.4 52.5 52.6 52.7 52.8 52.9 53.1 53.2 53.3
[277] 53.4 53.5 53.6 53.7 54.0 54.1 54.4 54.5 54.6 54.7 54.8 55.1
[289] 55.4 55.6 56.3 56.4 56.5 56.6 57.1 57.3 59.0 65.3 349.2 406.8
[301] 412.8 414.0 424.8 429.6 430.8 458.4 469.2 470.4 481.2 496.8 498.0 500.4
[313] 501.6 502.8 508.8 511.2 519.6 523.2 525.6 535.2 536.4 546.0 579.6 619.2
[325] 624.0
datos_preproc <-
datos |>
mutate(MedianAge = ifelse(MedianAge > 130, median(MedianAge), MedianAge))
unique(datos_preproc$MedianAge) |> sort() [1] 22.3 23.2 23.3 23.5 23.9 24.2 24.4 24.6 24.7 25.0 25.2 25.9 26.1 26.2 26.3
[16] 26.4 26.5 26.6 26.7 26.8 27.0 27.2 27.4 27.5 27.6 27.8 27.9 28.0 28.1 28.3
[31] 28.4 28.5 28.6 28.7 28.8 28.9 29.1 29.2 29.3 29.4 29.5 29.7 29.8 29.9 30.0
[46] 30.1 30.2 30.3 30.4 30.5 30.6 30.7 30.8 30.9 31.0 31.1 31.2 31.3 31.4 31.5
[61] 31.6 31.7 31.8 31.9 32.0 32.1 32.2 32.3 32.4 32.5 32.6 32.7 32.8 32.9 33.0
[76] 33.1 33.2 33.3 33.4 33.5 33.6 33.7 33.8 33.9 34.0 34.1 34.2 34.3 34.4 34.5
[91] 34.6 34.7 34.8 34.9 35.0 35.1 35.2 35.3 35.4 35.5 35.6 35.7 35.8 35.9 36.0
[106] 36.1 36.2 36.3 36.4 36.5 36.6 36.7 36.8 36.9 37.0 37.1 37.2 37.3 37.4 37.5
[121] 37.6 37.7 37.8 37.9 38.0 38.1 38.2 38.3 38.4 38.5 38.6 38.7 38.8 38.9 39.0
[136] 39.1 39.2 39.3 39.4 39.5 39.6 39.7 39.8 39.9 40.0 40.1 40.2 40.3 40.4 40.5
[151] 40.6 40.7 40.8 40.9 41.0 41.1 41.2 41.3 41.4 41.5 41.6 41.7 41.8 41.9 42.0
[166] 42.1 42.2 42.3 42.4 42.5 42.6 42.7 42.8 42.9 43.0 43.1 43.2 43.3 43.4 43.5
[181] 43.6 43.7 43.8 43.9 44.0 44.1 44.2 44.3 44.4 44.5 44.6 44.7 44.8 44.9 45.0
[196] 45.1 45.2 45.3 45.4 45.5 45.6 45.7 45.8 45.9 46.0 46.1 46.2 46.3 46.4 46.5
[211] 46.6 46.7 46.8 46.9 47.0 47.1 47.2 47.3 47.4 47.5 47.6 47.7 47.8 47.9 48.0
[226] 48.1 48.2 48.3 48.4 48.5 48.6 48.7 48.8 48.9 49.0 49.1 49.2 49.3 49.4 49.5
[241] 49.6 49.7 49.8 49.9 50.0 50.1 50.2 50.3 50.4 50.5 50.6 50.7 50.8 50.9 51.0
[256] 51.1 51.2 51.3 51.4 51.5 51.7 51.8 51.9 52.0 52.1 52.2 52.3 52.4 52.5 52.6
[271] 52.7 52.8 52.9 53.1 53.2 53.3 53.4 53.5 53.6 53.7 54.0 54.1 54.4 54.5 54.6
[286] 54.7 54.8 55.1 55.4 55.6 56.3 56.4 56.5 56.6 57.1 57.3 59.0 65.3
- `deathRate`: entre 0 e infinito.
- `medianIncome`: entre 0 e infinito.
- `popEst2015`: entre 0 y 8.000.000.000 (habitantes del mundo xd)
- `povertyPercent`: entre 0 y 100 (al ser porcentaje)
`studyPerCap`: entre 0 e infinito
`MedianAge`: entre 0 y 130 años (por poner un tope)
`AvgHouseholdSize`: entre 0 e infinito (si existe alguna finca enorme)
- `PercentMarried`: entre 0 y 100 (al ser porcentaje)
- `PctHS25_Over`: entre 0 y 100 (al ser porcentaje)
- `PctUnemployed16_Over`: entre 0 y 100 (al ser porcentaje)
- `ctPrivateCoverage`: entre 0 y 100 (al ser porcentaje)
- `BirthRate`: entre 0 e infinito.
a la hora de sustituir usa lo mediana porque es menos sensible a valores extremos como un 500 y pico que había
Ejercicio 2 (0.75 puntos)
En caso de encontrar alguna variable cualitativa, procesa los datos como consideras para incluir su información de manera que pueda ser usada en la futura regresión. Recuerda: si tenemos una variable cualitativa que toma k valores, debemos crear k-1 nuevas variables numéricas (de un tipo muy concreto) en su lugar.
# region cualitativa
unique(datos$region)[1] "West" "South" "Midwest" "Northeast"
# hay 4 valores posibles: creamos 3 variables nuevas binarias
datos_preproc <-
datos_preproc |>
mutate('West' = ifelse(region == "West", 1, 0),
'South' = ifelse(region == "South", 1, 0),
'Midwest' = ifelse(region == "Midwest", 1, 0)) |>
select(-region)region es cualitativa con 4 posibles calores: creamso 3 nuevas variables y quitamos region
Ejercicio 3 (1.5 puntos)
Ejecuta el código que consideres para realizar una selección inicial de variables de manera que se mitiguen los efectos adversos de la colinealidad. Comenta todo lo que consideres sobre. Tras ello separa las observaciones en dos datasets distintos: un dataset train y un segundo dataset test. Para ello realiza un muestreo aleatorio de manera que train contenga el 80% de los datos y test el 20% restante. IMPORTANTE: no cambies la semilla.
#correlacion con deathrate
#datos_preproc |> correlate() |> focus(deathRate)
#datos_preproc |> cor() |> corrplot(method='ellipse')
#quitamos predictoras muy pococ correladas para simplificarnos los calculos
datos_2<- datos_preproc |> dplyr::select(-c(studyPerCap,MedianAge,AvgHouseholdSize))
set.seed(12345)
split<-initial_split(datos_2, prop=0.8 )
train<-training(split)
test<-testing(split)por tomar una medida quitamos las variables poco correladas con la objetivo deathrate, en este caso las menores al |0.1| aunque tambien podriamos considerar quitar popEst2015
Ejercicio 4 (1 punto)
No usaremos el dataset test hasta el final para evaluar las predicciones (con un dataset que el modelo no conoce). Con el dataset train ejecuta el ajuste de regresión lineal multivariante saturado (habiendo hecho lo pedido en los ejercicios anteriores) y comenta de manera MUY RESUMIDA la salida (SOLO lo que puedas interpretar hasta ahora, y al grano, que se te acaba el tiempo)
ajuste_saturado<-lm(train, formula=deathRate~.)
ajuste_saturado |> summary()
Call:
lm(formula = deathRate ~ ., data = train)
Residuals:
Min 1Q Median 3Q Max
-108.344 -11.959 0.627 11.555 158.574
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.243e+02 1.210e+01 10.280 < 2e-16 ***
medIncome -3.574e-05 7.566e-05 -0.472 0.636724
popEst2015 -9.133e-07 1.350e-06 -0.676 0.498794
povertyPercent 6.144e-01 1.812e-01 3.391 0.000708 ***
PercentMarried -3.148e-01 9.774e-02 -3.221 0.001296 **
PctHS25_Over 1.267e+00 8.404e-02 15.079 < 2e-16 ***
PctUnemployed16_Over 1.551e+00 1.863e-01 8.324 < 2e-16 ***
PctPrivateCoverage 6.610e-02 8.496e-02 0.778 0.436627
BirthRate -5.372e-01 2.338e-01 -2.297 0.021683 *
West -8.133e+00 2.195e+00 -3.704 0.000217 ***
South 1.038e+01 1.953e+00 5.317 1.15e-07 ***
Midwest 3.136e+00 1.956e+00 1.603 0.109030
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.94 on 2425 degrees of freedom
Multiple R-squared: 0.3765, Adjusted R-squared: 0.3737
F-statistic: 133.1 on 11 and 2425 DF, p-value: < 2.2e-16
# comprobamos colinealidad entre predictoras con VIF
check_collinearity(ajuste_saturado)# Check for Multicollinearity
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
medIncome 4.33 [4.04, 4.64] 2.08 0.23 [0.22, 0.25]
popEst2015 1.19 [1.15, 1.26] 1.09 0.84 [0.80, 0.87]
PercentMarried 2.28 [2.15, 2.43] 1.51 0.44 [0.41, 0.47]
PctHS25_Over 1.75 [1.66, 1.86] 1.32 0.57 [0.54, 0.60]
PctUnemployed16_Over 2.00 [1.89, 2.12] 1.41 0.50 [0.47, 0.53]
PctPrivateCoverage 4.17 [3.89, 4.47] 2.04 0.24 [0.22, 0.26]
BirthRate 1.07 [1.04, 1.14] 1.04 0.93 [0.88, 0.96]
West 3.07 [2.87, 3.28] 1.75 0.33 [0.31, 0.35]
South 4.78 [4.46, 5.13] 2.19 0.21 [0.19, 0.22]
Midwest 4.29 [4.01, 4.61] 2.07 0.23 [0.22, 0.25]
Moderate Correlation
Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
povertyPercent 6.78 [6.30, 7.30] 2.60 0.15 [0.14, 0.16]
train2<- train|> select(-povertyPercent)
ajuste_saturado2<-lm(train2, formula=deathRate~.)
ajuste_saturado2 |> summary()
Call:
lm(formula = deathRate ~ ., data = train2)
Residuals:
Min 1Q Median 3Q Max
-103.288 -11.893 0.604 11.863 158.564
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.569e+02 7.375e+00 21.275 < 2e-16 ***
medIncome -1.690e-04 6.480e-05 -2.608 0.009175 **
popEst2015 -1.060e-06 1.352e-06 -0.784 0.433138
PercentMarried -4.777e-01 8.531e-02 -5.599 2.40e-08 ***
PctHS25_Over 1.228e+00 8.342e-02 14.721 < 2e-16 ***
PctUnemployed16_Over 1.669e+00 1.834e-01 9.096 < 2e-16 ***
PctPrivateCoverage -5.346e-02 7.747e-02 -0.690 0.490157
BirthRate -5.073e-01 2.342e-01 -2.167 0.030357 *
West -8.090e+00 2.200e+00 -3.677 0.000241 ***
South 1.117e+01 1.943e+00 5.745 1.03e-08 ***
Midwest 3.588e+00 1.956e+00 1.835 0.066667 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.99 on 2426 degrees of freedom
Multiple R-squared: 0.3736, Adjusted R-squared: 0.371
F-statistic: 144.7 on 10 and 2426 DF, p-value: < 2.2e-16
check_collinearity(ajuste_saturado2)# Check for Multicollinearity
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
medIncome 3.16 [2.96, 3.38] 1.78 0.32 [0.30, 0.34]
popEst2015 1.19 [1.14, 1.25] 1.09 0.84 [0.80, 0.87]
PercentMarried 1.73 [1.64, 1.83] 1.32 0.58 [0.55, 0.61]
PctHS25_Over 1.72 [1.63, 1.82] 1.31 0.58 [0.55, 0.61]
PctUnemployed16_Over 1.93 [1.82, 2.05] 1.39 0.52 [0.49, 0.55]
PctPrivateCoverage 3.45 [3.23, 3.69] 1.86 0.29 [0.27, 0.31]
BirthRate 1.07 [1.04, 1.13] 1.03 0.93 [0.88, 0.96]
West 3.07 [2.87, 3.28] 1.75 0.33 [0.31, 0.35]
South 4.71 [4.39, 5.06] 2.17 0.21 [0.20, 0.23]
Midwest 4.27 [3.99, 4.59] 2.07 0.23 [0.22, 0.25]
consideramos un VIF mayor a 5 demasiado por lo que eliminamos povertyPercent, recomprobamos y parece que las variables predictoras no estan significativamente correladas entre si
Ejercicio 5 (1.75 puntos)
Realiza una selección de modelos en base a los criterios BIC y AIC. Recuerda que para evitar incompatibilidad con otros paquetes, no debes hAcer library(MASS) sino MASS::… Comenta e interpreta todo lo que sepas de las salidas generadas. ¿Cuál penaliza más? ¿Por qué? ¿Qué ventajas tiene la selección de BIC? Interpreta los coeficientes de ambos modelos y comenta las diferencias (en caso de que existiesen)
#BIC
MASS::stepAIC(ajuste_saturado2, direction='both',k=log(nrow(train2)))Start: AIC=15137.43
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + PctPrivateCoverage + BirthRate + West +
South + Midwest
Df Sum of Sq RSS AIC
- PctPrivateCoverage 1 230 1172938 15130
- popEst2015 1 297 1173005 15130
- Midwest 1 1627 1174335 15133
- BirthRate 1 2269 1174977 15134
- medIncome 1 3287 1175995 15136
<none> 1172708 15137
- West 1 6535 1179243 15143
- PercentMarried 1 15155 1187863 15161
- South 1 15955 1188663 15163
- PctUnemployed16_Over 1 39994 1212702 15211
- PctHS25_Over 1 104760 1277467 15338
Step: AIC=15130.11
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + West + South + Midwest
Df Sum of Sq RSS AIC
- popEst2015 1 251 1173189 15123
- Midwest 1 1537 1174475 15126
- BirthRate 1 2125 1175064 15127
<none> 1172938 15130
- West 1 6310 1179248 15135
- medIncome 1 6658 1179596 15136
+ PctPrivateCoverage 1 230 1172708 15137
- PercentMarried 1 15328 1188266 15154
- South 1 16986 1189924 15157
- PctUnemployed16_Over 1 48377 1221315 15221
- PctHS25_Over 1 104776 1277714 15331
Step: AIC=15122.84
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
BirthRate + West + South + Midwest
Df Sum of Sq RSS AIC
- Midwest 1 1623 1174812 15118
- BirthRate 1 2108 1175297 15119
<none> 1173189 15123
- West 1 6201 1179390 15128
- medIncome 1 7272 1180460 15130
+ popEst2015 1 251 1172938 15130
+ PctPrivateCoverage 1 184 1173005 15130
- PercentMarried 1 15080 1188269 15146
- South 1 17441 1190630 15151
- PctUnemployed16_Over 1 48126 1221315 15213
- PctHS25_Over 1 109109 1282298 15332
Step: AIC=15118.41
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
BirthRate + West + South
Df Sum of Sq RSS AIC
- BirthRate 1 1742 1176554 15114
<none> 1174812 15118
+ Midwest 1 1623 1173189 15123
+ popEst2015 1 337 1174475 15126
+ PctPrivateCoverage 1 98 1174714 15126
- medIncome 1 9521 1184333 15130
- PercentMarried 1 13708 1188520 15139
- West 1 25366 1200178 15163
- South 1 31422 1206234 15175
- PctUnemployed16_Over 1 47178 1221990 15207
- PctHS25_Over 1 107638 1282450 15324
Step: AIC=15114.22
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
West + South
Df Sum of Sq RSS AIC
<none> 1176554 15114
+ BirthRate 1 1742 1174812 15118
+ Midwest 1 1257 1175297 15119
+ popEst2015 1 308 1176246 15121
+ PctPrivateCoverage 1 26 1176528 15122
- medIncome 1 8811 1185365 15125
- PercentMarried 1 15444 1191998 15138
- West 1 25515 1202069 15159
- South 1 32512 1209067 15173
- PctUnemployed16_Over 1 47741 1224296 15203
- PctHS25_Over 1 109486 1286040 15323
Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + West + South, data = train2)
Coefficients:
(Intercept) medIncome PercentMarried
1.541e+02 -2.141e-04 -4.624e-01
PctHS25_Over PctUnemployed16_Over West
1.233e+00 1.694e+00 -1.072e+01
South
8.706e+00
# salida ultima iteracion:lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + West + South, data = train2)
ajuste_BIC <-lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + West + South,
data = train2)
ajuste_BIC |> summary()
Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + West + South, data = train2)
Residuals:
Min 1Q Median 3Q Max
-105.425 -12.091 0.501 11.766 159.447
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.541e+02 5.944e+00 25.933 < 2e-16 ***
medIncome -2.141e-04 5.020e-05 -4.266 2.07e-05 ***
PercentMarried -4.624e-01 8.187e-02 -5.648 1.82e-08 ***
PctHS25_Over 1.233e+00 8.202e-02 15.038 < 2e-16 ***
PctUnemployed16_Over 1.694e+00 1.706e-01 9.930 < 2e-16 ***
West -1.072e+01 1.476e+00 -7.259 5.21e-13 ***
South 8.706e+00 1.062e+00 8.194 4.02e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22 on 2430 degrees of freedom
Multiple R-squared: 0.3715, Adjusted R-squared: 0.37
F-statistic: 239.4 on 6 and 2430 DF, p-value: < 2.2e-16
#AIC
MASS::stepAIC(ajuste_saturado2, direction='both',k=2)Start: AIC=15073.65
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + PctPrivateCoverage + BirthRate + West +
South + Midwest
Df Sum of Sq RSS AIC
- PctPrivateCoverage 1 230 1172938 15072
- popEst2015 1 297 1173005 15072
<none> 1172708 15074
- Midwest 1 1627 1174335 15075
- BirthRate 1 2269 1174977 15076
- medIncome 1 3287 1175995 15078
- West 1 6535 1179243 15085
- PercentMarried 1 15155 1187863 15103
- South 1 15955 1188663 15105
- PctUnemployed16_Over 1 39994 1212702 15153
- PctHS25_Over 1 104760 1277467 15280
Step: AIC=15072.13
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + West + South + Midwest
Df Sum of Sq RSS AIC
- popEst2015 1 251 1173189 15071
<none> 1172938 15072
- Midwest 1 1537 1174475 15073
+ PctPrivateCoverage 1 230 1172708 15074
- BirthRate 1 2125 1175064 15074
- West 1 6310 1179248 15083
- medIncome 1 6658 1179596 15084
- PercentMarried 1 15328 1188266 15102
- South 1 16986 1189924 15105
- PctUnemployed16_Over 1 48377 1221315 15169
- PctHS25_Over 1 104776 1277714 15279
Step: AIC=15070.65
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
BirthRate + West + South + Midwest
Df Sum of Sq RSS AIC
<none> 1173189 15071
- Midwest 1 1623 1174812 15072
+ popEst2015 1 251 1172938 15072
+ PctPrivateCoverage 1 184 1173005 15072
- BirthRate 1 2108 1175297 15073
- West 1 6201 1179390 15082
- medIncome 1 7272 1180460 15084
- PercentMarried 1 15080 1188269 15100
- South 1 17441 1190630 15105
- PctUnemployed16_Over 1 48126 1221315 15167
- PctHS25_Over 1 109109 1282298 15285
Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + West + South + Midwest,
data = train2)
Coefficients:
(Intercept) medIncome PercentMarried
1.536e+02 -2.012e-04 -4.713e-01
PctHS25_Over PctUnemployed16_Over BirthRate
1.238e+00 1.705e+00 -4.851e-01
West South Midwest
-7.770e+00 1.149e+01 3.564e+00
# salida ultima iteracion:lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + BirthRate + West + South + Midwest, data = train2)
ajuste_AIC <- lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + BirthRate + West + South + Midwest,
data = train2)
ajuste_AIC |> summary()
Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + West + South + Midwest,
data = train2)
Residuals:
Min 1Q Median 3Q Max
-103.450 -11.926 0.524 11.870 158.525
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.536e+02 6.226e+00 24.673 < 2e-16 ***
medIncome -2.012e-04 5.187e-05 -3.879 0.000108 ***
PercentMarried -4.713e-01 8.436e-02 -5.587 2.58e-08 ***
PctHS25_Over 1.238e+00 8.236e-02 15.027 < 2e-16 ***
PctUnemployed16_Over 1.705e+00 1.709e-01 9.980 < 2e-16 ***
BirthRate -4.851e-01 2.323e-01 -2.089 0.036829 *
West -7.770e+00 2.169e+00 -3.582 0.000347 ***
South 1.149e+01 1.912e+00 6.008 2.16e-09 ***
Midwest 3.564e+00 1.945e+00 1.833 0.066953 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.98 on 2428 degrees of freedom
Multiple R-squared: 0.3733, Adjusted R-squared: 0.3712
F-statistic: 180.8 on 8 and 2428 DF, p-value: < 2.2e-16
BIC:
deathRate = 1.541e+02+-2.141e-04*medIncome + -4.624e-01*PercentMarried + 1.233e+00*PctHS25_Over +1.694e+00* PctUnemployed16_Over + -1.072e+01*West + 8.706e+00*South
AIC:
deathRate = 1.536e+02 +
medIncome * -2.012e-04 +
PercentMarried * -4.713e-01 +
PctHS25_Over * 1.238e+00 +
PctUnemployed16_Over * 1.705e+00 +
BirthRate * -4.851e-01 +
West * -7.770e+00 +
South * 1.149e+01 +
Midwest * 3.564e+00 +
Como vemos enl BIC penaliza más, quita dos predictoras mas que AIC, las predictoras birthrate y midwest , que no son buenas: cambia precision por senicillez.
Ejercicio 6 (1.5 puntos)
Compara la calidad de los 3 modelos (saturado, BIC y AIC) en train. Quédate con uno (justifica) y comprueba si se cumplen las hipotesis necesarias (fase de diagnosis SOLO en ese elegido). ¿Se cumplen las hipótesis? Argumenta porqué, más allá de lo que valga luego la bondad de ajuste, es importante que se cumplan.
compare_performance(ajuste_saturado2, ajuste_AIC,ajuste_BIC)# Comparison of Model Performance Indices
Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | R2 | R2 (adj.) | RMSE | Sigma
--------------------------------------------------------------------------------------------------------------------
ajuste_saturado2 | lm | 21991.6 (0.154) | 21991.7 (0.151) | 22061.1 (<.001) | 0.374 | 0.371 | 21.936 | 21.986
ajuste_AIC | lm | 21988.6 (0.690) | 21988.6 (0.691) | 22046.5 (0.013) | 0.373 | 0.371 | 21.941 | 21.982
ajuste_BIC | lm | 21991.5 (0.156) | 21991.6 (0.158) | 22037.9 (0.987) | 0.372 | 0.370 | 21.972 | 22.004
Obviamente el AIC ofrece menor AIC y el BIC menor BIC y el saturado presenta el mayor AIC y BIC.
Como es de esperar tambien el saturado ofrece mayor R2 (infromación eplicada) que AIC y BIC y AIC mayor R2 que BIC.
El BIC penaliza más, con menos variables pero el AIC ofrece mas informacion (apenas) a mas parametros aunque estos sean malos y por lo general es preferible elegir un modelo mas simple, especialmente la diferencia de informacion explicada es tan pequeña (0.002 ó 0.001 en R2 adj).
–>BIC es el mejor en este caso (sencillo a similar precisión)
Es necesario realizar inferencia estadística delos estimadores de los parámetros para comprobar que no solo funcionan en esta muestra, si no en otra y se pueda aplicar a toda la población
test linealidad -> OK
check_model(ajuste_BIC)chisq.test(ajuste_BIC$residuals,ajuste_BIC$fitted.values, simulate.p.value = TRUE, B=250)Pearson's Chi-squared test with simulated p-value (based on 250 replicates) data: ajuste_BIC$residuals and ajuste_BIC$fitted.values X-squared = 5936532, df = NA, p-value = 0.003984test homocedasticidad -> NO
necesitamos que la varianza del error sea finita y constante, tal que no varíe según aumenta o disminuye la X.
performance::check_heteroscedasticity(ajuste_BIC)Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).test normalidad de erroes -> NO
ε∼N(0,σ2r)
performance::check_normality(ajuste_BIC)Warning: Non-normality of residuals detected (p < .001).
independencia de errores -> OK
performance::check_autocorrelation(ajuste_BIC)OK: Residuals appear to be independent and not autocorrelated (p = 0.622).
Ejercicio 7 (0.75 puntos)
Asume que se cumplen las hipótesis (en caso de que no se cumpliesen) e interpreta todo lo que sepas sobre la segunda tabla de la salida de la regresión (la de inferencia de los parámetros). Elimina variables si su efecto no fuese significativo.
ajuste_saturado2 |> summary()
Call:
lm(formula = deathRate ~ ., data = train2)
Residuals:
Min 1Q Median 3Q Max
-103.288 -11.893 0.604 11.863 158.564
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.569e+02 7.375e+00 21.275 < 2e-16 ***
medIncome -1.690e-04 6.480e-05 -2.608 0.009175 **
popEst2015 -1.060e-06 1.352e-06 -0.784 0.433138
PercentMarried -4.777e-01 8.531e-02 -5.599 2.40e-08 ***
PctHS25_Over 1.228e+00 8.342e-02 14.721 < 2e-16 ***
PctUnemployed16_Over 1.669e+00 1.834e-01 9.096 < 2e-16 ***
PctPrivateCoverage -5.346e-02 7.747e-02 -0.690 0.490157
BirthRate -5.073e-01 2.342e-01 -2.167 0.030357 *
West -8.090e+00 2.200e+00 -3.677 0.000241 ***
South 1.117e+01 1.943e+00 5.745 1.03e-08 ***
Midwest 3.588e+00 1.956e+00 1.835 0.066667 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.99 on 2426 degrees of freedom
Multiple R-squared: 0.3736, Adjusted R-squared: 0.371
F-statistic: 144.7 on 10 and 2426 DF, p-value: < 2.2e-16
Ahora que se ‘cumplen las hipótesis’ podemos interpretar el p-valor, R2 y otros estadisticos
Ejercicio 8 (1.5 puntos)
Por último, utiliza el dataset test para, aplicando los 3 modelos entrenados (saturado, AIC y BIC), predecir las observaciones de test (con cada uno de los modelos. Con las observaciones de test construye un nuevo dataset de TRES columnas: la variable objetivo, la variable predictora, la predicción y el modelo usado (saturado, BIC, AIC). Calcula el \(R^2\) en el dataset de test para cada modelo. ¿Cuál funciona mejor en test? ¿Qué % de información explica cada modelo? ¿Cuánto vale su varianza residual?
pred_test_S <-
tibble("objetivo" = test$deathRate,
"prediccion" = predict(ajuste_saturado2, newdata = test),
"split" = "test")
pred_test_A <-
tibble("objetivo" = test$deathRate,
"prediccion" = predict(ajuste_AIC, newdata = test),
"split" = "test")
pred_test_B <-
tibble("objetivo" = test$deathRate,
"prediccion" = predict(ajuste_BIC, newdata = test),
"split" = "test")
pred_train_S <-
tibble("objetivo" = train$deathRate,
"prediccion" = ajuste_saturado2$fitted.values,
"split" = "train")
pred_train_A <-
tibble("objetivo" = train$deathRate,
"prediccion" = ajuste_AIC$fitted.values,
"split" = "train")
pred_train_B <-
tibble("objetivo" = train$deathRate,
"prediccion" = ajuste_BIC$fitted.values,
"split" = "train")
prediccion <- bind_rows(pred_test_S, pred_test_A, pred_test_B, pred_train_S, pred_train_A, pred_train_B)
prediccion# A tibble: 9,141 × 3
objetivo prediccion split
<dbl> <dbl> <chr>
1 144. 155. test
2 237. 206. test
3 207. 191. test
4 210. 208. test
5 156. 183. test
6 140. 193. test
7 169. 199. test
8 158. 159. test
9 176. 185. test
10 187. 173. test
# ℹ 9,131 more rows